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Abstract 



In this paper, we address the problem of safety verification of interval hybrid systems 
in which the coefficients are intervals instead of explicit numbers. A hybrid symbolic- 
^ i numeric method, based on SOS relaxation and interval arithmetic certification, is pro- 

]^ I posed to generate exact inequality invariants for safety verification of interval hybrid 

systems. As an application, an approach is provided to verify safety properties of non- 
1/^ . polynomial hybrid systems. Experiments on the benchmark hybrid systems are given 

^SJ ' to illustrate the efficiency of our method. 

o ■ 

cn 

1. Introduction 

^ ! As a tool of modelling cyber-physical systems, hybrid systems are dynamical systems gov- 
I erned by interacting discrete and continuous dynamics. The continuous dynamics of a hybrid 
system is specified by differential equations, and for discrete transitions, the hybrid system 
changes state instantaneously and possibly discontinuously. Among the most important re- 
search issues in formal analysis of hybrid systems are safety, i.e., deciding whether a given 
property holds in all the reachable states, and its dual problem reachability, i.e., deciding if 
there exists a trajectory starting from the initial set that reaches a state satisfying the given 
property. Due to the infinite number of possible states in state spaces, safety verification 
and reachability analysis of hybrid systems presents a challenge. For general (exact) hybrid 
systems, some well-established techniques [26, 7, 16, 21, 32, 22, 35, 34] based on invariant 
generation have been proposed for safety verification of the systems. However, when apply- 
ing these techniques, one can not avoid numerical errors or may suffer from high complexity. 
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To take advantage of the efficiency of numerical computation and the error-free property 
of symbohc computation, we proposed in [36] a hybrid symbohc-numeric method via exact 
sums-of-squares (SOS) representation to construct differential invariants for continuous dy- 
namic systems, and generalized in [15, 37] the idea for safety verification of polynomial hybrid 
systems. 

A common assumption made on hybrid systems is that the coefficients of the involved 
equations are specific values. In practice, however, due to the increasing complexity of 
modern systems, some disturbance and modeling errors may be contained in the system 
description, and, in addition, there may be noisy and inexact data involved in the realistic 
problem. All these factors may contribute to inexactness of the data used to describe the 
hybrid systems. To take this uncertainty into account, it would be more reasonable and 
appropriate to use intervals rather than concrete but inexact data to represent the hybrid 
systems. This motivates us to introduce the notion of interval polynomial hybrid systems, by 
which we mean the differential equations in hybrid systems are represented as polynomials 
with interval coefficients. 

In this paper, we consider safety verification of interval polynomial hybrid systems, i.e., 
deciding whether none of trajectories of an interval hybrid system starting from the initial 
set can enter some unsafe regions in the state spaces. In [37] we applied a symbolic-numeric 
computation method, based on bilinear matrix inequality (BMI) solving and exact SOS 
polynomials representations, to deal with exact safety verification for polynomial hybrid 
systems. In this paper, we extend the techniques in [37] to generate exact invariants for 
verifying interval hybrid systems. The idea lies in applying interval arithmetic to verify 
positive semidefiniteness of interval matrices and existence of solutions to interval polynomial 
equations. As an apphcation, we apply the above approach to verify safety of non-polynomial 
hybrid systems by relaxing continuous dynamics of non-polynomial forms to those of interval 
polynomial forms, and then studying safety of the latter system whose set of trajectories 
contains that of the original non-polynomial system. 

The contributions of our paper are as follows. First, an approach is proposed to verify 
safety property of an interval hybrid system, therefore, safety property is guaranteed for an 
arbitrary hybrid system within the given interval system. Moreover, our approach can gener- 
ate exact invariants instead of approximate ones, overcoming the unsoundness of verification 
caused by numerical errors [20]. And in comparison with some symbolic approaches based 
on qualifier elimination, our approach is more efficient and practical, because parametric 
polynomial optimization problem based on SOS relaxation can be solved in polynomial time 
theoretically. Second, a key problem we consider in safety verification is that of determining 
nonnegativity of interval multivariate polynomials, which is a fundamental problem in real 
algebraic geometry. Thirdly, for a non-polynomial function, we propose a rigorous polynomial 
approximation method to compute its approximate polynomial with polynomial lower and 
upper bounds of the interpolation error. Compared with the classical Taylor approximation, 
the polynomial bounds we give is much sharper. 

The rest of the paper is organized as follows. In Section 2, we introduce some notions 
related to interval hybrid systems. Section 3 is devoted to determining nonnegativity of 
interval multivariate polynomials. In Section 4, two techniques which combine SOS relaxation 
with interval arithmetic are proposed to generate invariants of interval hybrid systems with 
small and large radii, respectively. As an application, safety verification of non-polynomial 
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hybrid systems is discussed in Section 5. Section 6 concludes the paper. 

2. Interval Hybrid Systems and Safety Verification 

Let us first review some notions of general hybrid systems [9, 32]. 

Definition 1 (Hybrid System). A hybrid system is a tuple H : (y, L, T, 0, 2^, ^o) with 

• V = {xi, ...,Xn}, a set of real-valued system variables; 

• L, a finite set of locations; 

• if) & L, the initial location; 

• T , a set of transitions. Each transition r : (£, Qr, Pr) G T consists of a ^relocation i G 
L, a postlocation i' G L, the guard condition over V , and an assertion pr over V U 
V representing the next-state relation, where V = {x[, ...,x'j^} denotes the next-state 
variables; 

• Q, an assertion specifying the initial condition; 

• T>, a map that associates each location I E L to a differential rule (a.k.a. a vector 
fieldj 'D{t), an autonomous system Xi = fe,i(y) for each Xi G V, written briefly as 
X = f£(x); 

• , a map that sends I E L to a location invariant '^{t), an assertion over V . 

In reality, due to measuring errors or disturbance, the data involved in the systems may 
be inaccurate. It is then reasonable to consider hybrid systems in which some data are given 
as interval estimates rather than specific values, the so-called interval hybrid systems. Similar 
to Definition 1, an interval hybrid system IH is defined to be a tuple 

(F,L,r,e,[p],vi/,4), 

where V , L, T, 6, Eq are the same as in Definition 1, while [V] represents a map sending 
each location £ G L to an interval differential rule [^^(^)] of the form 

Xi = [/]£,i(x) i = l,...,n, 

by which we mean [/]^^j(x) is a real function with interval coefficients; for brevity, we write 
[!'(£)] as X = [f]£(x); For more details on interval arithmetic, please refer to Appendix A. 

A hybrid system H : (V, L, T, 6, V, \E^, Iq) is said to be within an interval hybrid system 
IH : (V, L, T, 6, [D], io) if /£,i(x) G [f]£,i{x) for each i e L and z = 1, . . . , n, or written 
briefly as G 

In this paper, we will mainly study safety verification of interval hybrid systems. Recall 
that a hybrid system is said to be safe if none of the trajectories starting from any state in 
the initial set can evolve to an unsafe region. Similarly, given a prespecified unsafe region 
Xu C M*^, an interval system IH : (V,L,T, 0, ['D],\E',£o) is said to be safe if every hybrid 
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system within IH is safe. This is to say, none of the trajectories of interval hybrid system 
IH starting from any state in the initial set can evolve to Xu, or, equivalently, any state in 
Xu is not reachable. 

Recall that an invariant of a hybrid system H is an over- approximation of all the reachable 
states of the system H. Since generating invariants of arbitrary form for hybrid systems is 
computationally hard, the usual technique is to compute inductive invariants. It is shown 
in [37] that safety verification of general hybrid systems can be reduced to finding inductive 
invariants (a.k.a. barrier certificates in [22]) of hybrid systems, as described in the following 
theorem. 

Theorem 1. [[22], [36]] Let H : (V, L, T, 0,P, \&,^o) be a general hybrid system. Suppose 
that for each location i ^ L, there exists a function (pe{'x.) such that 

(i) 0hv^^o(x)>O; 

(ii) </?f(x) > A g{i,i') A p{i,i') \= v^f(x') > 0, for any transition {£,£',g,p) going out of i; 

(iii) (pe{'x.) > A \E'(£) |= </3^(x) > 0, where ^^(x) denotes the Lie-derivative of ipiix.) along 
the vector field V{i), i.e., ^^(x) = Yh=i ^feA^)- 

Then (fei'^) > is an (inductive) invariant of the hybrid system H at location L If, moreover, 

(iv) X„(£) h Vti?^) < for any i e L, 
then the safety of the system H is guaranteed. 

The notion of inductive invariants can be generalized for interval hybrid systems, as 
defined in the following 

Definition 2 (Inductive invariant). For an interval hybrid system IH : {V, L, T, O, [T>], Eq) , 
an inductive assertion map X of IH is a map that associates with each location i ^ L an 
assertion X{i) that holds initially and is preserved by all discrete transitions and continuous 
flows o/IH. More formally, the map X satisfies the following requirements: 

[Initial] e h^(4)- 

[Discrete Consecution] For each discrete transition r : {£,£', gr, Pr) starting from a state 
satisfying X{i), taking r leads to a state satisfying X{i'), i.e., X{i) A gr A pr \= X{i') 
where X{i') represents the assertion X{i) with the current state variables Xi 's replaced 
by the next state variables x\ 's, respectively. 

[Continuous Consecution] For location i ^ L and states xi), X2) such that X2 
evolves from Xi according to any differential rule T>{i) G ['P](^), if^i \= X{i) then X2 |= 
X{i). 

The difference between inductive invariants of interval hybrid systems and those of general 
hybrid systems lies in that for continuous consecution, any differential rule contained in the 
interval differential rule must be considered. Then Theorem 1 can be modified for verifying 
safety of interval hybrid systems, as described in the following. 
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Theorem 2. LetlH : (V, L, T, 0, [P],\l/,^o) an interval hybrid system. Suppose that for 
each i ^ L, there exists a function (feix.) satisfying the conditions (i-ii) in Theorem 1, and 

(iii') <^^(x) > 0A^(£) 1= y3^(x) > 0, here <~pi{^ denotes the Lie-derivative of (pe{x.) along any 
differential rule !){£) E [!){£)], i.e., Lpe{x) = X;r=i ^kii^)' M a^V feA^) ^ IfiAi^)- 

Then feix.) > is an (inductive) invariant of the interval hybrid system IH at location i. 
If, moreover, the condition (iv) in Theorem 1 is satisfied, then the safety of the system IH 
is guaranteed. 

In our preceding papers [15, 37], a symbolic-numeric method based on SOS relaxation, 
Gauss-Newton refinement and rational vector recovery techniques is proposed to generate 
polynomial inequality invariants </?£(x) > at each location £ E L for general polynomial 
hybrid systems. This method can not be applied directly on interval hybrid systems. In 
the sequel, we will combine BMI solving with interval arithmetic to compute polynomial 
invariants </?£(x) > which satisfy conditions in Theorem 2. For brevity, we will abuse the 
notation (feix) to represent both the polynomial (fe{'^) and the invariant V9f(x) > 0. 

3. Nonnegativity of Interval Polynomials 

To determine whether a polynomial inequality v^£(x) > is an invariant of an interval hybrid 
system, by Theorem 2 (iii') it suffices to decide whether a multivariate polynomial V9^(x) with 
interval coefficients is positive semidefinite. In the sequel, we will call a polynomial with 
interval coefficients an interval polynomial. Denote by IM[x] the set of interval multivariate 
polynomials in x. The first problem to be investigated is the following 

Problem 1. Given an interval polynomial [ip]{'x.) G I]R[x], verify whether it is positive 
semidefinite, or the validity of the interval inequality 

[V-llx) > 0, Vx G M". 

It is well known that the problem of testing positive semidefiniteness of real polynomials is 
NP-hard (when the degree is at least four). As stated in Appendix B, a sufficient condition 
for a multivariate polynomial to be positive semidefinite is that there exists an SOS polyno- 
mial (or rational function) representation. In [10, 11, 19], some symbolic-numeric methods 
were proposed to determine whether a multivariate polynomial ipi^) with rational coeffi- 
cients is positive semidefinite by computing its exact SOS representations, or equivalently, to 
determine if there exists a symmetric matrix W G M^'^'^ satisfying exactly 

ip{x) = m(x)^ ■ W ■ m(x) and W hO, (1) 

where W ^ denotes that W is positive semidefinite. These methods cannot be applied 
directly to verifying positive semidefiniteness of an interval polynomial ['?/'](x) G IM[x], since 
there are infinitely many polynomials in the interval, and it is impossible to provide certifi- 
cates of SOS representations for infinitely many polynomials in [V'](x). For Problem 1, we 
will only prove existence of SOS representations for polynomials in [^/'](x). This problem can 
be further distinguished into two cases according to the radii of the coefficient intervals: the 
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coefficient intervals of [^](x) > are all smaller (resp. larger) than the given threshold. In 
the sequel, we will describe how to deal with the former case, and the latter case will be 
discussed in subsection 4.1. 

Let [W] be an interval matrix such that [W] >z 0, i.e., every matrix within [W] is positive 
semidefinite. If for any polynomial ipix.) within [■?/'] (x), there exists a matrix W G [W] 
such that the condition (1) holds exactly, then we have [^/'](x) > 0. Thus the first case of 
Problem 1 can be transformed into the problem of finding an interval matrix [W] >z for 
M(x). 

Suppose that there exists an approximate SOS decomposition of the mid-point function 
midV^(x) e 

m.idip{x.) ^ m(x)^ ■ W ■ m(x) (2) 

where W ^ 0. Having W, we will consider how to compute an interval matrix [W] ^ of 
minimal radius, such that [W] contains W and for any iplx) G ['?A](x) there always exists a 
matrix W G [W] satisfying the condition (1) exactly. Considering whether the matrix W is 
of full rank, there are two cases to be addressed. 

3.1. I? is of full rank 

Suppose that W in (2) is of full rank numerically, namely, the minimal eigenvalue of W is 
greater than the given tolerance r > 0. Let 

[W] ■.= W+[AW] 

be an interval matrix perturbed from W where [AW] G IM^^'^'. If, for any 'i/'(x) G ['i/'](x), 
there exists a matrix AW G [AW] which satisfies 

7/;(x) = m(x)^-(l? + Aiy)-m(x), (3) 

and W + AW >: exactly, then we have [^/'](x) > 0. Since W is positive definite and of full 
rank, according to matrix perturbation theory we have W + [AW] ^ as long as the radius 
of interval matrix [AW] is small enough. 

We ffist consider how to construct an interval matrix [AW] with small radius, which 
satisfies the condition (3). Comparing the coefficients of terms on both sides of (3) gives rise 
to the following underdetermined linear system with the entries of AW as unknowns w: 

A-w = [v], 

where A G M*^'' with s G Z"*" and r = M^±li^ w G IW is a vector composed of columnwise 
entries of the symmetric matrix AW, and [v] G IW is the coefficient vector of the interval 
polynomial [^/'](x) — m(x)"^ ■ W ■ m(x). Our goal is to compute a minimal 2— norm interval 
vector w satisfying A ■ w = [v] . The above problem is then transformed into the following 
interval least squares problem: 

S = min{||w||2 : A ■ w = v for some v G [v]}. 
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Using the method [31] for solving interval linear systems, we can obtain a solution [w'] G IMJ" 
of S and therefore the associated solution [A 14^] of (3) of minimal radius. Then the remaining 
task is to verify whether the interval matrix 14^ + [AH^] is positive semidefinite. The following 
theorem provides such a computational criterion. 

Theorem 3. [27, Theorem 4] Let [W] be a symmetric interval matrix and [W] = [W — 
AW, W + AW] be its midpoint-radius form. Suppose that p{AW) is the spectral radius of 
AW and XminiW) is the minimum eigenvalue of W . If p{AW) < \min{W), then [W] is 
positive semidefinite. Moreover, if p{AW) < Amm(W^) then [W] is positive definite. 

We give an example to illustrate the above method. 
Example 1. Verify ['?/'] (x) > where 

[V^](x) = 0.9574 - 1.9362x1 - 0.3404x2 + [1.1852, 1.2593]x? 

- [0.4237, 0.4576]xi xa + [1.125, 1.2083]x^. 

For the mid-point function mid^/'(x), we compute its approximate Gram matrix representa- 
tion midip{'x) ^ m(x)'^ ■ W ■ m(x) where 

-0.1702 \ 
-0.2203 . 
1.1667 J 

It is easy to check that W is of full rank. By solving an associated interval linear system, we 
obtain the symmetric interval matrix [W] as follows: 

/ [0.9574,0.9575] [-0.9681,-0.9680] [-0.1703,-0.1702] \ 

[-0.9681,-0.9680] [1.1851,1.2593] [-0.2289,-0.2118] . 
y [-0.1703,-0.1702] [-0.2289,-0.2118] [1.1388,1.1945] y 

For the midpoint-radius form of [W], we obtain 0.0422 = p{AW) < Xmin{W) = 0.0461. 
According to Theorem 3, [W] is positive definitive, which proves [^/'](x) > 0. □ 

3.2. W is singular 

When the matrix W is singular or near to a singular matrix, the perturbed matrix of W may 
not be positive semidefinite. Therefore, the method in subsection 3.1 does not apply to the 
case where W is numerically singular. 

By expanding the quadratic representation, the equation (2) can be rewritten as 

1=1 ^ a 



m X 



/ 1 \ 

Xi 

U2 / 



/ 0.9574 -0.9681 
W = -0.9681 1.2222 
. -0.1702 -0.2203 
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where / is the rank of W. Next we will verify, for each tp^x) G [?/'](x), there exist qi^a £ ^ 
such that 



holds exactly. Let q be a vector composed of all the qi^a- Comparing the terms of both sides 
of (4) gives rise to a nonlinear system of the form 



where F : W ^ W with r the size of q, and [v] G IR^ is an interval vector consisting of 
coefficients in ['?/'] (x). Note that F{0) = 0. Hence, the problem of determining [^/'](x) > is 
equivalent to that of verifying existence of real roots of the underdetermined interval nonlinear 
system (5). The latter problem can be solved in two ways: one is based on existence of real 
roots for particular interval square nonlinear systems, and the other for particular interval 
underdetermined nonlinear systems. The details of these two methods are given in Appendix 



Remark 1. If we find a verified real solution to system (5), then 'i^(x) > for each iplx) G 
['?/'] (x). However, the opposite is not true, i.e., even if ['?/'](x) > it is not guaranteed that 
the above methods can prove existence of real roots of (5). 

4. Safety Verification of Interval Hybrid Systems 

In this section, we study how to verify safe properties of an interval hybrid system. Two 
techniques will be used depending on the radii of the occurred intervals in the given interval 
hybrid system. If the radii of the intervals are all larger than a given threshold, we trans- 
form the interval hybrid system into an uncertain hybrid system by replacing the intervals 
with some uncertainties and then generalize the method in [15, 37], which is based on SOS 
relaxation and rational vector recovery, to compute exact invariants of the uncertain hybrid 
system. If the radii of the involved intervals are all less than the given threshold, we will 
apply the interval verification approach in Section 3. For the more general case, when the 
interval hybrid system contains both intervals of radii smaller than and those of radii larger 
than the given threshold, the above two techniques will be combined. For simplification, we 
will only consider the two special cases respectively in subsections 4.1 and 4.2. 

4.1. Safety Verification of Interval Hybrid Systems With Large Radii 

Let IH : {V, L, T, O, [T>], i^) be an interval hybrid system. Suppose that the radii of the 
intervals in the interval differential rules [D] are all greater than a given threshold e, say 
e = 0.1. Then some new parameters ui,...,Ut will be introduced to replace the interval 
coefficients, to convert IH into an uncertain hybrid system Hu with u = {ui, . . . ,Ut), for 
which Theorem 1 can be extended to handle safety verification. 

Denote by [u] = [u, u] G IM* the interval coefficient vector composed of all the interval 
coefficients occurred in [D], where u = {ui, . . . ,mJ and u = {ui, . . . ,Ut)- To remove the 




(4) 



F(q) - [v] = 0, 



(5) 



C. 
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intervals [u] in IH, we introduce a vector u G of uncertainties with the constraints 

■di{u) = {ui -Ui)(ui -Ui) >0, i = l,...,t. 

For the uncertain hybrid system Hu, we predetermine a template </?(x) = Cq,x" of 
polynomial invariants with the given degree d, where x° = ■ ■ ■ a;"", a = {ai, . . . , a„) G 
Z>Q with Yyi=i — Cq G M are parameters. For each location £ G L, we write 

(pii'x.) = cj ■ T(x), where T(x) is the (column) vector of all terms in xi, . . . ,Xn with total 
degree < d, and G M'^, with = ("^'^), is the coefficient vector of ^eix). For clarity, 
we write </5^(x) as (/?f(x, c^). Similar to Theorem 1, the problem of computing the invariants 
(fiei'x.) of the uncertain hybrid system Hu can be translated into the following problem 

' find ce G W, E L 
s.t. e 1= v9f„(x, > 0, 

(^,(x,c,) >OA(?(£,f)Ap(£,f) h<^^'(x',QO >0, (6) 
M^, > A A 79(u) > h 0£(x, u, q) > 0, 

^ 1= (^^(x,q) < 0, 

where ^^(x, u, q) = J^i^i |^ ' fiA^^^)- Without loss of generality, we consider a simpler 
form of (6): 

find ceW 
s.t. V5i(x,c)>0, 

<^3(x,c) >0 |=(^2(x,c) >0, ^ ^ 

(/P5(X, u, c) > h V54(x, u, c) > 0, 

where the coefficients of the polynomials yjj's are affine in c, for i = 1, . . . , 5. By Appendix 
B, the problem (7) can be further transformed into the following polynomial parametric 
optimization problem 

' find c G W 

s.t. V9i(x, c) = mi(x)^ ■ ■ mi(x), 
< V52(x, c) = m2(x)^ ■ . ni2(x) + (m3(x)'^ ■ W^^^ ■ m3(x)) • V93(x, c), 

V94(x, u, c) = m4(x, u)^ ■ W ■ m4(x, u) + (m5(x, u)^ ■ W^^^ ■ m5(x, u)) ■ v^slx, u, c), 
VrW^O, z = l,...,5, 

(8) 

which involves both LMI and BMI constraints. As stated in [37], a Matlab package PENBMI 
solver [13], which combines the (exterior) penalty and (interior) barrier method with the 
augmented Lagrangian method, can be applied directly on the BMI program, and alterna- 
tively, an iterative method can be applied by fixing W^^^ and c alternatively, which leads to 
a sequential convex LMI problem. 

Since the SDP solvers in Matlab is running in fixed precision, the above techniques yield 
numerical vector c and numerical positive semidefinite matrices W^^\ . . . , W^^\ which satisfy 
the constraints in (8) approximately. We will apply the symbolic-numeric method proposed in 
[37] to obtain exact solutions to (8). The idea is as follows. We first convert W^^'^ and W^^'^ to 
the nearby rational positive semidefinite matrices W^^^ and W^^\ respectively, by nonnegative 
truncated PLDL^P"'"-decomposition, in which all the diagonal entries of the corresponding 
diagonal matrix are preserved to be nonnegative. Then, using modified Newton refinement 
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and rational vector recovery techniques, we can recover the rational vector c and the ra- 
tional positive semidefinite matrices W^^\ W^'^\ W^^'^ from the numerical c, H^W, VTt^', VTt^', 
respectively, such that the constraints in (8) hold exactly. For more details, please refer to 
[37]. 

4.2. Safety Verification of Interval Hybrid systems with Small Radii 

In this subsection, we will consider interval hybrid systems with small radii interval coef- 
ficients, namely, the radii of the involved intervals are all smaller than the given threshold 
e. For such interval hybrid systems, the method described in subsection 4.1 via introduc- 
ing uncertainties may suffer from high complexity especially when solving the parametric 
optimization problem (8). Instead, we will consider how to generate invariants of IH by de- 
termining nonnegativity of interval polynomials: we first compute candidate invariants with 
rational coefficients, then employ the interval computation method presented in Section 3 to 
certify that the candidate invariants satisfy the conditions in Theorem 2 exactly. 

Suppose that [-D](^) of IH is given by x = [f£](x) for i E L. Choosing the midpoints of 
the interval coefficients of [f£](x) yields a mid-point vector midf£(x) and an associated general 
hybrid system H with the vector field x = midff(x), for £ E L. Then the symbolic- numeric 
technique in [37] can be used to generate invariants of H as follows. Let us predetermine a 
polynomial template (fe{'x.) > of invariants of H with deg(pe{x.) = d. By Theorem 1, the 
problem of computing V9^(x) can be translated into the following problem 

' find Q G G L 

s.t.e h v5£o(x,Qo) > 0' 

< ^,(x,q) > OA^?(£,f) Ap(£,f) h<^^'(x',c,0 > 0, (9) 

V5^(x, q) > A h midv9^(x, q) > 0, 
^ X„(£) 1= y?£(x,Q) < 0, 

where mid(/3^(x, q) = X]r=i ' By use of BMI solving and modified Newton 

refinement, we can obtain the refined numerical solutions to (9). With the refined vector 
for i E L, we then apply rational vector recovery technique to obtain a polynomial (/^^(x, c^) 
with rational coefficients. Clearly, V9^(x, c^) can be seen as a candidate invariant of the interval 
hybrid system IH. 

In the following, we will determine whether V9^(x, c^) satisfies the conditions of invariants 
of interval hybrid system IH in Theorem 2 exactly, i.e., 

e \= v?^o(x,Cfj > 0, 

(^,(x, c,) > A g{i, f ) A p{i, £') 1= <^,,(x', ce>) > 0, 

<^,(x, c,) > A h [Vi] (x, ce) > 0, ^ ^ 

X^{i) \= <^^(x,c^) < 0, 

where [v9f](x, q) = Y17=i |^ ' [fe,i]i^) interval polynomial. Observing in (10), all the 

constraints except the third one are exact constraints. And the SOS-based method presented 
in subsection 4.1 can be used to determine satisfiability of the exact constraints. To handle 
the third constraint in (10), we now consider how to determine satisfiability of polynomial 
inequalities with interval coefficients. More generally, we consider the following problem 

V'i(x,c) >0 1= [^/;2](x,c) >0, (11) 
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where ['?/'2](x, c) G I]R[x]. Let mid^/'(x, c) G [^/'2](x, c) be the mid-point function of [?/'2](x, c). 
Then BMI solver and modified Gauss-Newton refinement can yield the numerical positive 
semidefinite matrices W^^'^ and W^'^\ which satisfy the following condition approximately 

midV'lx, c) ^ m2(x)^ ■ l^l^] . ^^^(x) + (mi(x)^ ■ ■ mi(x)) ■ ipi{^, c). (12) 

Converting W^^'^ to a nearby rational positive semidefinite matrix W^^'^ by nonnegative trun- 
cated PLDL"'"P"'"-decomposition, the condition (12) becomes 

midV'(x, c) - (mi(x)^ ■ W^W ■ mi(x)) ■ ipi{^, c) ^ m2(x)'^ ■ W^^^ ■ m2(x). (13) 

Let [1^2] (x, c) be an interval polynomial such that 

[V^2](x,a) = M(x,e) - (mi(x)^ ■ . ni^(x)) . ^i(x,c). 

bmce [^1 h 0, it suffices to prove satisfiability of (11) when [^/'2](x, c) is nonnegative. Remark 
that (13) is an approximate SOS decomposition of ['?/'2](x, c). The nonnegativity of [ijj2] can 
be verified by computing the corresponding interval matrix [W2], either using the method 
in subsection 3.1 if W^t^' is of full rank, or by proving existences of real roots of the interval 
nonlinear system, as explained in subsection 3.2. 



4.3. Experiments 







. *2 _ 





In the following, some examples will be given to illustrate our method on safety verification 
of interval hybrid systems. 

Example 2. Consider the classical two-dimensional system given in [12, 22], whose coeffi- 
cients are approximated and described by the following intervals 

[0.99, 1.01]x2 
■ [0.96, 1.04]xi + [0.32, 0.347]a;? - [0.98, 1.02]x2 

We will verify that all trajectories of the system starting from the initial set 

e = {(xi, X2) G : (xi - 1.5)^ + xl< 0.25} 

will never enter the unsafe region 

Xu = {{xi, X2) G : (xi + If + (X2 + 1)' < 0.16}. 

Set the threshold e = 0.1 Clearly, all the radii of involved intervals are less than this 
threshold. Applying the method in subsection 4.2, we obtain the following verified invariant 
with rational coefficients 



y?(x) 



151 152 



Xi 



99 99 

which guarantees the safety of the original system. 



62 106 

— X2 H 

33 99 



X1X2 + -x^, 
y 



□ 
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Example 3. Consider a Moore- Greitzer model of a jet engine with stabilizing feedback op- 
erating in the no-stall mode [2]. In this model, the origin is translated to a desired no-stall 
equilibrium. The dynamic system takes the following form: 

xi = [-1.1, -0.9]a;2 - fx^ - 
X2 = [2.98,3.02]xi -X2. 

The problem is to verify that all trajectories of the system starting from the initial set 

e = {(xi, xa) e : (xi - 1)^ + xj < 0.04} 

will never reach the unsafe set 

= {(xi, xa) G : (xi + 1.8)^ + xl < 0.16}. 

Set the threshold e of radii to be 0.1. Then a new uncertainty u is introduced to replace 
the interval [—1.1,-0.9]. Combine the methods in subsections 4.1 and 4.2 to deal with 
the uncertain interval system, and we obtain the following verified invariant with rational 
coefficients 



(p{Xi,X2) 



2231 652 



274 



46 



10 



328 123 

which guarantees the safety of the original system. 



Xo X, H XiXo 

123 41 ^ 41 



1649 
"98? 



25 



□ 



Example 4. Figure 4 gives a predator-prey hybrid system [24] with interval coefficients: 

A(X) = /2(X) = 



-Xl + [0.99, 1.01]xiX2 
[0.875, 1.2]x2 - X1X2 



Suppose the system starts in location ii with an initial state in 

e = {(xi, X2) e : (xi - 0.8)2 ^ _ Q < Q.Ol}. 
We want to verify that the system never reach the states in 

g(l,2) : (.T2 - 0.875)(Z2 - 0.9) < 




p{2, 1) : (2:1 - 0.7)^ + {X2 - O.T f < 0.01 



Figure 1: Hybrid system of Example 4 



X^{ii) = {(a;i, X2) e : 0.8 < Xi < 0.9 A 0.8 < X2 < 0.9}. 
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Set the threshold e of radii to be 0.1. Then a new uncertainty u is introduced to replace 
the interval [0.875, 1.2]. Applying the above method on the resulting uncertain interval hybrid 
system IH^, we obtain the following verified invariants with rational coefficients 

411 346 397 49 ^ 

(/JifXi, Xo) = \ Xi H Xo Xn, 

^ ^ ' 995 995 995 199 ^ 



, 556 151 986 22 , 

CZJolXi , Xo) = Xi Xo H Xn, 



for locations l\ and £2, respectively, which ensures the safety of the original hybrid system. 
□ 



5. Safety Verification of Non-polynomial Hybrid system 

As an application of the method in Section 4 for safety verification for interval hybrid systems, 
we will consider how to verify safety of non-polynomial hybrid systems. 

Let H : {V, L,T,Q,T>,'$ , Eq) be a hybrid system where the initial condition O, location 
invariants "^{i), the guard condition and reset relation in each transition r G T are semi- 
algebraic sets, whereas the continuous systems in the differential rules 'D{i), contain some 
non-polynomial terms in x. For such a non-polynomial hybrid system H, we will first trans- 
form it into an uncertain interval hybrid system IH through polynomial approximation on 
the non-polynomial terms, such that H is within IH. This implies that the safety of IH 
suffices to guarantee the safety of H, and then the method in Section 4 can be applied to the 
former problem. 

Assume that the location invariant "^{i) is a compact set for each location i. Consider 
the continuous dynamics of a hybrid system H at location i: 

s 

Xi = /i(x) = /io(x) + ^ /ij(x)0ij(x), z = 1, . . . , n, (14) 

where x takes values in C M"-, /jj(x) are polynomials for j = 0, 1, . . . , s, and 0ij(x) are 
non-polynomials for j = 1, . . . , s. We will approximate the functions 0ij(x) with polynomials 
Qiji'x) for i = 1, . . . ,n and j = 1, . . . , s. Let /ijj be the bound of |0ij(x) —gij(x.)\ for all x G 
namely, 

I (x) - Qij (x) I < flij , for all X G ^ (£) . (15) 

Making use of the relation (15) for each location i E L, we can construct an interval poly- 
nomial hybrid system IH : {V,L,T,Q, [D], where the interval differential rule ['P(^)] 
given by 

s 

Xi = [/i](x) = /io(x) + fij{:s^){gij{x) + [-/iij, /^ij]), i = 1, . . . , n, (16) 

j=i 

enclosures the non-polynomial system (14) in H, that is, /i(x) G [/i](x) for all x G \l/(£). 

The key point of the above idea is to compute an approximate polynomial and the associ- 
ated bound for the given non-polynomial function. For a non-polynomial function 0(x) with 
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X G we will compute the approximate polynomial g{x.) G M[x] with a verified bound 

/i G M+, such that 

|0(x)-(7(x)|</i,Va;Gvi>(£), 

and the bound fi is as small as possible. 

A classic method of polynomial approximation is Taylor expansion. In this paper, to ob- 
tain a tighter error bound, multivariate polynomial interpolation [6] is applied to compute an 
approximate polynomial with the error bound. Furthermore, the technique of oversampling is 
explored to get better approximate polynomials, i.e., the number of the interpolation points 
is greater than that of the terms of the target polynomial 5'(x). Given the interpolation 
points, the approximate polynomial g{x.) can be obtained by solving a least squares problem. 
Specifically, predetermine a polynomial template of g{x.) with a given degree r: 

(7(x)=c^-r(x), (17) 

where T(x) is the (column) vector consisting of all terms in Xi, . . . , x„ with total degree < r, 
and c G M*^, with u = ("^^), is the coefficient vector of 5'(x). We then construct a mesh M 
on "^{i) with a small spacing s G IR+, and compute Uj = 0(vj) G M for 1 < j < m at mesh 
points {vi, V2, Vm}. Let the coefficient vector c of 5f(x) be unknowns. We can construct a 
linear system 

A-c = y, (18) 

where A = (T(vi)"^, T(v2)"^, T(vm)"^)^ is of size m x z/ with m> u. By solving the above 
overdetermined system, we obtain (/(x, c) as the approximation of (/)(x) with x G '^{^)- Having 
5f(x, c), one will compute the verified error bound namely, |0(x) — 5f(x, c)| < ;U, Vx G '^{t). 

Lemma 1. [38, Theorem 3] Let K (Z he a convex polyhedron, and Vi,V2, ...Vm and d 
he the vertices and diameter of K respectively. Suppose : K ^ M. is a continuous and 
differential function on K , then for all ai, a2, ...a^ G such that ai + a2 + ... + = 1, we. 
have 

Tl 

\ip{x) - (ai^(Vi) + a2ip{V2) + ... + a ))l < Pd, 

n + 1 

where (3 = sup^gx || \/ ip{'x)\\. 

For the error function r(x) = 0(x) — g{'x,c), we will estimate its bound with x in the 
mesh M by the following theorem. 

Theorem 4. Suppose that s and {vi, V2, v^} are the mesh spacing and mesh points of 
M, respectively. Lei /xq = niax{r(vi), r(v2), r(vm)}, and f3' = sup^^J^J \\ \/ r(x.)\\, then for 
all x G M, 

71 

\r{^)\<^j(3's + iio. 
n + 1 

Proof. We know that r(x) is a continuous and differential function on M. Thus, according 
to Lemma 1, for all Oil , fl2 , . . .flm G such that ai + ^2 + ... + = 1, 

71 

|r(x) - (air(vi) + a2r(v2) + ... + amr(v^))| < 
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IT 



Figure 2: Approximate e^, — 2 < x < 2 by g{x,c) + [— /x, /i] (solid line: e^, dot line: g{x,c) ± fi). 



Then, we have 



ft Tl 

k(x)| < — —rf^'s + |(air(vi) + a2r{\2) + ••• + a^r{\„^))\ < — —rf^'s + fio- 
n + 1 n + 1 



□ 



Example 5. Consider the function = with \1/ : —2 < a; < 2. We want to compute a 
polynomial g{x) and the associated verified error bound fi such that 



g{x)\ < fx, -2<x <2. 



First, we construct a mesh M on \E' with the spacing s = |. For a polynomial of the form 
g{x, c) = Co + Ci X + C2 + C3 x^, it is easy to find an approximate polynomial 

^(x, c) = 0.9173 + 0.9562X + 0.6797x^ + 0.2117x1 

According to Theorem 4, we can also compute the error bound fi = 0.2937. The results are 
as shown in Figure 2. □ 

Stated as above, once we obtain an interval polynomial hybrid system IH from H through 
polynomial approximation such that H is within IH, the method in Section 4 can be used 
to verify safety of IH, which ensures safety of H. The following example is presented to 
illustrate our method for safety verification of a non-polynomial hybrid system. 



Example 6. Consider the following two-tanks hybrid system [25] depicted in Figure 3 with 
/i(x) = 



1 - ^/xT 
xi - Jx^ 



/2(X) 



1 - a/Xi - X2 + 1 



^Xx - X2 + 1 - a/x^ 



where x = (xi,X2) denotes the liquid levels. In [25], the authors verified that the system 
starting in location with an initial state in 

e = {(xi,X2) G : (xi - 5.5)2 ^ _ q 25)2 < 0.0625} 

15 




Figure 3: Hybrid system of Example 6 

never reach the states of 

= {(xi, X2) E : (xi - 4.25)2 + (xa - 0.25)^ < 0.0625}. 
Here, we enlarge both radii of initial and unsafe regions to 0.49, that is, 
e = {(xi, X2) G : (xi - 5.5)2 ^ (^^^ _ q 25)2 < o.2401} 

and 

X„(£i) = {(xi, X2) e : (xi - 4.25)2 + (xs - 0.25)^ < 0.2401}, 

and consider again safety verification of the given system. We first compute an interval 
polynomial system given by [/i](x) and [/2](x) to enclosure the original system where 

[0.1658, 0.173] - 0.3377x1 + 0.0114x? 
[0.6465, 0.8615] + 0.3377xi - 1.7115x2 - 0.0114x2 + 0.8241x2 



[/i](x 
and 

[/2](X)^ 



-[0.1204, 0.132] - 0.3316x1 + 0.3319x2 + 0.0135x2 - 0.0269xiX2 + 0.0137x| 
[0.6716, 0.6898] + 0.3316xi - 0.9576x2 - 0.0135x? + 0.0269xiX2 + 0.0572xi 

We obtain the following verified invariants with rational coefficients 



^ , , 1069 145 367 121 2 160 242 2 

09i(x) = Xi X2 H Xi H X1X2 H Xo, 

^ ^ ' 994 142 497 497 ^ 497 497 2' 

^ , , 9621 20 899 989 , 349 1487 , 

= W - 71^^ + 497^^ + 994^^ + 497^^^^ ' W^2, 
which satisfy the conditions in Theorem 2 exactly. Therefore, the safety of the original hybrid 
system is verified. □ 

The above approach can be easily extended to the case of uncertain non-polynomial hybrid 
systems, by which we mean the continuous dynamics at each location ^ are given by uncertain 
non-polynomial systems of the form 

s 

X, = /,(x, 9) = /,o(x, e) + Y, /u(x, ^)0^,(x), for 1 < z < n, (19) 

i=i 

where 6' G $ C M* is a vector of uncertainty. The following example demonstrates how to 
apply the above approach to verify safety of an uncertain non-polynomial system. 
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Example 7. Consider an uncertain non-polynomial system given in [5]: 



xi = -xi + X2 + |(ei - 1), 

X2 = —Xi — X2 + 6X1X2 + Xi COSXi, 

for —2 < Xi,X2 < 2 and 0.98 <9< 1.2. This system starts with an initial state in 

e = {(xi, X2) eM.^ ■.x\ + xl< 0.25}, 

and we want to verify that the system never reach the states of 

Xu = {{xu X2) G : + 1.5)^ + (xs + 1.5)^ < 0.16}. 

To prove the safety of this non-polynomial system, we first compute interval polynomials to 
approximate the non-polynomial terms e^^ and cos Xi occurred in this system. Based on the 
above techniques, we obtain the following interval polynomial system 

= [-0.1882, 0.1055] - 0.5219xi + X2 + 0.33985x? + 0.10585x?, 
X2 = [-0.2067, 0.0875]xi - X2 + ^XiXa - 0.3594x?, 

which enclosures the original system. For the above interval hybrid system, we obtain the 
following verified invariant with rational coefficients 

, 343 31 25 49 . 17 55 n 

C5 X = 1 Xi H Xo Xi Xi Xo Xn, 

' 32 6 48 32 ^ 48 32 
which guarantees the safety of the original system. 



6. Conclusion 

In this paper, a hybrid symbolic-numeric method, based on SOS relaxation and interval 
arithmetic certification, is proposed to generate exact inequality invariants for safety veri- 
fication of interval hybrid systems. As an application, one approach is provided to verify 
safety property of non-polynomial hybrid systems. More precisely, we apply a rigorous poly- 
nomial approximation method to compute an interval polynomial system, which contains the 
non-polynomial system, and then compute the exact invariant of the corresponding interval 
polynomial system to verify the safety property of the original system. Experiments on the 
benchmark hybrid systems illustrate the efficiency of our algorithm. 



References 

[1] Alefeld, G., and Herzberger, J. Introduction to interval computations. 1983. 

[2] Aylward, E., Parrilo, P., and Slotine, J. Stability and robustness analysis 
of nonlinear systems via contraction metrics and sos programming. Automatica 44i 8 
(2008), 2163-2170. 



17 



[3] BOCHNAK, J., COSTE, M., AND RoY, M. Real Algebraic Geometry, vol. 36. Springer 
Verlag, 1998. 

[4] Chen, X., and Womersley, R. Existence of solutions to systems of underdetermined 
equations and spherical designs. SI AM Journal on Numerical Analysis 44i 6 (2006), 
2326-2341. 

[5] Chesi, G. Estimating the domain of attraction for non-polynomial systems via LMI 
optimizations. Automatica 45, 6 (2009), 1536-1541. 

[6] Gasca, M., and Sauer, T. On the history of multivariate polynomial interpolation. 

Journal of computational and applied mathematics 122, 1 (2000), 23-35. 

[7] GULWANI, S., and Tiwari, a. Constraint-based approach for analysis of hybrid 
systems. In CAV (2008), vol. 5123 of LNCS, Springer, pp. 190-203. 

[8] Hansen, E. Bounding the solution of interval linear equations. SI AM journal on 
numerical analysis 29, 5 (1992), 1493-1503. 

[9] Henzinger, T. The theory of hybrid automata. In Proceedings of the 11th Annual IEEE 
Symposium on Logic in Computer Science (1996), IEEE Computer Society, pp. 278-292. 

[10] Kaltofen, E., Li, B., Yang, Z., and Zhi, L. Exact certification of global op- 
timality of approximate factorizations via rationalizing sums-of-squares with floating 
point scalars. In Proceedings of the International Symposium on Symbolic and Algebraic 
Computation (New York, NY, USA, 2008), ISSAC, ACM, pp. 155-164. 

[11] Kaltofen, E., Li, B., Yang, Z., and Zhi, L. Exact certification in global polynomial 
optimization via sums-of-squares of rational functions with rational coefficients. Journal 
of Symbolic Computation 4^ (2012), 1-15. 

[12] Khalil, H. Nonlinear Systems, 2rd ed. New Jewsey, Prentice hall, 1996. 

[13] KOCVARA, M., and Stingl, M. PENBMI user's guide (version 2.0). Available at 
http://www.penopt.coni, 2005. 

[14] Krawczyk, R. Newton-algorithmen zur bestimmung von nullstellen mit fehler- 
schranken. Computing 4, 3 (1969), 187-201. 

[15] Lin, W., Wu, M., Yang, Z., and Zeng, Z. Exact safety verification of hybrid 
systems using sums-of-squares representation. Accepted for publication in SCIENCE 
CHINA Information Sciences, 15 pages, 2012. 

[16] Liu, J., Zhan, N., and Zhao, H. Computing semi-algebraic invariants for polyno- 
mial dynamical systems. In Proceedings of the International Conference on Embedded 
Software (EMS OFT) (2011), ACM, pp. 97-106. 

[17] LOFBERG, J. YALMIP: A toolbox for modeling and optimization in MAT- 
LAB. In Proceedings of the CACSD (Taipei, Taiwan, 2004). Available at 
http : // control . ee . ethz . ch/~ joloef /yalmip .php. 



18 



Parrilo, p. Structured Semidefinite Programs and Semialgebraic Geometry Methods 
in Robustness and Optimization. PhD thesis, Cahfornia Institute of Technology, 2000. 

Peyrl, H., and Parrilo, P. Computing sum of squares decompositions with rational 
coefficients. Theoretical Computer Science 409, 2 (2008), 269-281. 

Platzer, a., and Clarke, E. M. The image computation problem in hybrid systems 
model checking. In Hybrid Systems: Computation and Control, HSCC (2007), Springer, 
pp. 473-486. 

Platzer, A., and Clarke, E. M. Computing differential invariants of hybrid systems 
as fixedpoints. Form. Methods Syst. Des. 35, 1 (2009), 98-120. 

Prajna, S., Jadbabaie, A., and Pappas, G. A framework for worst-case and 
stochastic safety verification using barrier certificates. IEEE Transactions on Automatic 
Control 52, 8 (2007), 1415-1429. 

Prajna, S., Papachristodoulou, A., and Parrilo, P. SOSTOOLS: 
Sum of squares optimization toolbox for MATLAB, 2002. Available at 
http : //www. cds . caltech.edu/sostools. 

Ratschan, S., and She, Z. Safety verification of hybrid systems by constraint propa- 
gation based abstraction refinement. Hybrid Systems: Computation and Control (2005), 
573-589. 

Ratschan, S., and She, Z. Safety verification of hybrid systems by constraint 
propagation-based abstraction refinement. ACM Transactions in Embedded Comput- 
ing Systems 6, 1 (2007), 573-589. 

Rodriguez-Carbonell, E., and Tiwari, a. Generating polynomial invariants for 
hybrid systems. In Hybrid Systems: Computation and Control, HSCC (2005), vol. 3414 
of LNCS, pp. 590-605. 

Rohn, J. Positive definiteness and stability of interval matrices. SI AM Journal on 
Matrix Analysis and Applications 15, 1 (1994), 175-184. 

Rohn, J., and Kreinovich, V. Computing exact componentwise bounds on solutions 
of lineary systems with interval data is NP-hard. SIAM Journal on Matrix Analysis and 
Applications 16, 2 (1995), 415-420. 

Rump, S. On the solution of interval linear systems. Computing J^l, 3 (1992), 337-353. 

RUMP, S. Intlab-interval laboratory. Developments in Reliable Computing (1999). 

Rump, S. Verification methods: Rigorous results using fioating-point arithmetic. Acta 
Numerica (2010), 287-449. 

[32] Sankaranarayanan, S., Sipma, H., and Manna, Z. Constructing invariants for 
hybrid systems. Formal Methods in System Design 32 (2008), 25-55. 



19 



[33] Sturm, J. F. Using SeDuMi 1.02, a MATLAB toolbox for optimization over symmetric 
cones. Optimization Methods and Software 11/12 (1999), 625-653. 

[34] Sturm, T., and Tiwari, A. Verification and synthesis using real quantifier elimina- 
tion. In Proceedings of the International Symposium on Symbolic and Algebraic Com- 
putation, ISSAC (2011), ACM Press, pp. 329-336. 

[35] Tiwari, A. Approximate reachability for linear systems. In Hybrid Systems: Compu- 
tation and Control, HSCC (2003), vol. 2623 of LNCS, pp. 514-525. 

[36] Wu, M., AND Yang, Z. Generating invariants of hybrid systems via sums-of-squares 
of polynomials with rational coefficients. In Proc. 2011 Internat. Workshop on Symbolic- 
Numeric Comput. (New York, N. Y., 2011), ACM Press, pp. 104-111. 

[37] Yang, Z., Wu, M., and Lin, W. Exact verification of hybrid systems based on 
bilinear sos representation. Submitted, 19 pages, 2012. 

[38] ZENG, Z., AND ZHANG, J. A mechanical proof to a geometric inequality of zi- 
rakzadeh through rectangular partition of polyhedra (in Chinese). Journal of Systems 
Science and Mathematical Sciences 30, 11 (2010), 1430-1458. 

Appendix 

A. Interval Arithmetic 

Interval arithmetic [1] has been designed for automatically tackling roundoff errors of numer- 
ical computations. In this subsection, some notions about interval arithmetic are presented. 

Denote the closed intervals by [a], [b], etc. By convention, the left and right endpoints of 
a closed interval [a] are represented by a and a, respectively, i.e., 

[a] = {x G M : a<x<a} 

with a = mi[a] and a = sup [a]. Any real number a can also be regarded as an interval by 
identifying a with the "point interval" [a] with a = a = a. Such point intervals are also called 
degenerate intervals. Let 

mid([a]) := -{a + a) and rad[a] := -{a — a) 

be the midpoint and radius of the closed interval [a], respectively. Clearly, an interval can 
also be represented by its midpoint and radius. The set of all intervals over M is denoted 
by IM. The arithmetic operations +,—,*, ^ can be extended from M to IM in the usual set 
theoretic sense, and the bounds of the resulting intervals can be computed from the bounds 
of the operands, see [1] for details. 

By IM" and IM"^^" we denote the sets of real n- dimensional vectors and m x n matrices 
over IM, respectively. Elements of IM" are called interval vectors and denoted by [a], [b] 
and etc, and elements of IM™^*^ are called interval matrices and denoted by [A], [B] and etc. 
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Remark that interval vectors (resp. interval matrices) are sets of vectors (resp. matrices). 
For interval vectors and matrices, the notions of midpoints and radius, and the arithmetic 
operations are defined componentwise. 

By an interval linear system, we mean a system of the form 

[A]x=[b], (20) 

where [A] [b] G M" and x column vector of n unknowns. The 

set 

S = {x G M" : Ax = b for some A E [A],h E [h]} 

is called the solution set of the interval system (20). Many efficient algorithms are available 
in [8, 28, 29, 31] for obtaining guaranteed inclusions [31] for the solution set S. 

Let f : M" — > M" be a continuously differentiable function. Replacing the real vector x 
by an intervector [x] G IM" we thus obtain an interval extension [f] of f . By the inclusion 
property of interval arithmetic, the range of f over an interval is contained in its interval 
extension, i.e. {f(x) : x G [x]} C [f]([x]). To determine existence of solutions to the 
nonlinear system f(x) = 0, we will use the Krawczyk operators [14] based on Browder fixed 
points, which is defined as follows. Assume that [x] G IM" is an interval set satisfying x G [x], 
and C E R"^'^. The Krawczyk operator is defined as follows 

K{±, [x], f) = X - Cf (x) + (/ - C[f' ]([x]))([x] - x). 

In practical computation, C is usually chosen to be near the inverse of the Jacobian f (x). 

Theorem 5. [31] Under the above assumptions, if 

K{±, [x],f)Czni([x]), 

where int{\x\) is the topological interior of [x], then there exists a unique x* G -ft'(x, [x],f) 
such that f(x*) = 0. 

INTLAB is a MATLAB toolbox [30], which consists of interval calculations, and interval 
arithmetic for vectors and matrices. Many interval operations in this paper are implemented 
in MATLAB that uses the INTLAB package supporting rigorous real interval standard func- 
tions and interval least squares problem. 

B. Sum of Squares Relaxation 

We give a brief review on SOS optimization. More details will be found in [18]. Recall that 
a sufficient condition for determining E ]R[x] to be positive semidefinite is that there 

exists an SOS decomposition of ?/'(x): 

s 

^(x) = ^/,2(x), with/,(x)GR[x], (21) 

i=l 

or, equivalently, ^/'(x) can be represented in the Gram matrix form 

?/^(x) = m(x)"^ ■ W ■ m(x). 
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where TV is a real symmetric and positive semidefinite matrix over R, and m(x) is a vector 
of all monomials in M[x] with degree < |degr(x). Therefore the SOS program (21) can be 
further converted into the following Semidefinite programming (SDP) problem 

inf Trace(iy) 

s. t. ^(x) = m(x)'^ ■ W ■ m(x) ) (22) 

where Trace(iy) acts as a dummy objective function that is commonly used in SDP for 
optimization problem with no objective functions. Many Matlab packages of SDP solvers, 
such as SOSTOOLS [23], YALMIP [17], and SeDuMi [33], are available to solve the SDP 
problem (22) efficiently. 

The SOS programs have many applications, for example, in determining the nonnegativity 
of a multivariate polynomial over a semialgebraic set. Consider the problem of verifying 
whether the implication 

m 

/\(p,(x) >0) =^ g(x) >0 (23) 

i=l 

holds, where Pi(x) G ]R[x] for 1 < z < m and g(x) G M[x]. According to Stengle's Positivstel- 
lensatz, Schmiidgen's Positivstellensatz or Putinar's Positivstellensatz [3], if there exist SOS 
polynomials cJi G ]R[x] for i = 0, ...,m, such that 



q{x.) = ao(x) + ^(Ti(x)pi(x) 



1=1 



then the assertion (23) holds. Therefore, the existence of SOS representations provides a 
sufficient condition for determining the nonnegativity of g(x) over {x G M*^ : /\^^pj(x) > 0}. 



C. Existence of Real Roots for Underdetermined Interval Nonlinear 
Systems 

Consider a nonlinear system 

F(q) - [v] = 0. (24) 

where F : M'' — )■ a continuously differentiable function with r > s, r = Dim(q), and 
[v] G IW. To determine the existence of solution to system (24), we present two methods 
as follows. The idea of the first method is to transform the underdetermined interval system 
into the corresponding interval square nonlinear system by fixing some variables as constants, 
and then generalize Theorem 5 to verify the existence of real roots for this square nonlinear 
system, while the second method is to generalize the method in [4] to solve the interval 
underdetermined system (5). 

Firstly, suppose that q is an approximate solution of (24). Here we assume that the 
Jacobian matrix F'{q) at q is of full row rank. Column pivoting Qi?-decomposition for -F'(q) 
is applied to choose an index set B = {ki,k2, ■■■,ks} such that F^{q) G is nonsingular, 
that is, 

F'(q) = Q[R\*] with P G M"^", Q, i? G M'^^ 
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where P is a permutation matrix, Q is orthogonal and R is upper triangular. The permutation 
P arises from a greedy strategy to obtain maximum diagonal elements in R. Then, the set 
B can be taken as those components which are permuted to the first s positions by P. Thus, 
q can be separated into two parts q = (q_B,qAf), where N = {1,2, ...r}/B. Similar to the 
partition of q, we have q = (q^, q^). By use of the evaluations qN = qj^, (24) becomes the 
following interval square system 

G{qB) - [v] = 0, (25) 

where G(qB) = F(qB, cin) — c, and [v] = c + [v], and c is the constant vector of F(qB, qAr), 
i.e., c = F(0,qAr). 

Observing in (25), the interval coefficients only occur in the constant vector [v], and 
G{qB) is a real function from M'^ to W, meaning that the Jacobian matrix of (25) is the same 
as one exact square system G{qB) — v = 0, where v is a vector chosen randomly. Taking 
advantage of this property, it is easy to generalize Theorem 5 to verify the existence of real 
roots for (25). 

Theorem 6. Consider the system (25). Let [q^] G JMJ^ he such that q^ G [qs], and C G 
W'. If 

K{qB, [qs], G - [v]) = q^ - C{G{qB) - [v]) + (/ - CG"([qB]))([qB] - qs) C int{\qB\), 

(26) 

then there is a unique root q^ of (25) in [q^] for each v G [v]. 

Proof. If (26) holds, then we have, for each v G [v], 

K{qB, [qsl G - v) = qB - C{G{qB) - v) + (/ - CG'([qB]))([qB] - qs) C int([qB]). 

According to Theorem 5, for v there exists a unique root q^ of (25) in [q^]. Hence, for each 
V G [v], there exists a unique root for the system (25) if (26) holds. □ 

Alternatively, we also can apply another method provided in [4], to determine the exis- 
tence of real roots for the underdetermined system (5) directly. The only difference is that 
we need deal with a special interval underdetermined system while they worked on an exact 
one. For the same reason as in the above discussion, it is easy to generalize their method in 
[4] to deal with our problem. 

Suppose the Jacobian F\q) is of full row rank. Following [4], we apply the column 
pivoting Q-R-decomposition to choose an index set S = {ki,k2, ...,ks} such that -Fe(q) G M''^" 
is nonsingular. Then, define the function if : R*" — )■ by 

r HB{q) = qB-F'j,{qr\F{q)-w), 
\ HN{q) = qN- «(qAf - qw), 

where N = {1, 2, ...r}/B and a G (0, 1) is a constant. Obviously, if q* G M*" is a fixed point of 
H, that is H{q*) = q*, then we have F{q*) — v = with q^ = q^r. Choose two nonnegative 
numbers Vi and r2, we define the convex set 

[q] = {q G W : ||qB - qB|| < ri, Hq^v - q7v|| < ra}. 

Now, we can use the following theorem to determine the existence of solutions to the system 
(5). 
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Theorem 7. Consider the system (5). Suppose the Jacobian -F'(q) has full row rank, and 
that 

\\FM - Km < - q|l for q G [q]. 

There is a solution q* of (5) in [q] for each v G [v] if 

max ||F]j(q)-^(F(q) - v)|| + ||F^(q)-i (^Aln + r^^, + max ||F;,(q) ||r2) < r^. 

ve[v] 2 qG[q] 

Remark that Theorem 6 is a special case of Theorem 7 by setting r2 = 0. Compared with 
Theorem 7, the condition (26) in Theorem 6 is easy to verify in practice. 
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